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A study of a model rod-like polyelectrolyte molecule immersed into a monovalent or divalent 
electrolyte is presented. Results for the local concentration profile, mean electrostatic potential, 
charge distribution function and (^—potential are obtained from hypernetted-chain/mean spherical 
approximation (HNC/MSA) theory and compared with molecular dynamics (MD) simulations. As 
a particular case, the parameters of the polyelectrolyte molecule are mapped to those of a DNA 
molecule. Both, HNC/MSA and MD, predict the occurrence of overcharging, which is not present 
in the Poisson-Boltzmann theory. Father an excellent qualitative, and in some cases quantitative, 
agreement between HNC/MSA and MD is found. Oscillations observed in the mean electrostatic po- 
tential, local concentration profiles, as well as the curvature of the (^-potential are discussed in terms 
of the observed overcharging effect. Particularly interesting results are a very non-monotonic behav- 
ior of the (^-potential, as a function of the rod charge density, and the overcharging by monovalent 
counterions. 
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INTRODUCTION 

"Polyelectrolytes are polymers bearing ionizable 
groups, which, in polar solvents, can dissociate into 
charged polymer chains (macroions) and small counteri- 
ons" jlj. The combination of macromolecular properties 
and long-range electrostatic interactions results in an im- 
pressive variety of phenomena. It makes these systems 
interesting from a fundamental as well as a technologi- 
cal point of view. A thorough understanding of polyelec- 
trolytes has become increasingly important in biochem- 
istry and molecular biology. This is due to the fact that 
virtually all proteins, as well as the DNA, are polyelec- 
trolytes. Their interactions with each other and with the 
charged cell-membrane are still far from being fully un- 
derstood, which is partly due to the intricate coupling 
between ion distribution and chain conformation. 

A first approach to the problem is to fix the conforma- 
tion of the chain and to focus on a detailed description 
of the counterion distribution. Usually polyelectrolytes 
stretch due to the electrostatic repulsion of their charged 
groups. Moreover, many important polyelectrolytes have 
a large intrinsic stiffness (e.g., DNA, actin filaments or 
microtubules). Therefore, a rod-like conformation is an 
obvious first choice. The remaining problem of charged 
rods immersed into solution is much easier, but is still far 
from being exactly solvable. 

An additional approximation, which is frequently used 
in theoretical descriptions, is to completely integrate out 
the counterionic degrees of freedom. On a linearized 



mean-field level this yields a Debye-Hiickel-like theory 
characterized by a screened Coulomb potential between 
charged monomers. To obtain the correct physical prop- 
erties, one uses an effective Yukawa potential, which in 
turn requires adjustable parameters such as an effective 
polyelectrolyte radius and charge. In this way the depen- 
dence of the ionic structure, hereafter called electrical 
double layer (EDL), around two or more polyelectrolytes, 
as a function of the polyelectrolyte-polyelectrolyte dis- 
tance, is not accounted for. This information, however, 
is most relevant for the understanding of polyelectrolyte 
aggregation or self-assembling. This effect is more impor- 
tant for polyelectrolytes in low concentration added salt 
solution, since for this case the ionic screening is weak. 
This has been shown to be relevant for charged plates 
and charged spherical macroions |^ Q . 

A common further approximation assumes that the in- 
vestigation of a small sub-volume containing only one 
rod and its counterions will suffice to unveil much of the 
interesting physics. The main justification for this ap- 
proach is that the sub- volume has zero net charge. More- 
over, the counterions will also efficiently screen higher or- 
der multipoles. Hence, the interactions between two such 
sub- volumes, which are neglected when focusing on just 
one rod, will be fairly weak. This approximation is called 
cylindrical cell model and it provides the framework for 
our simulation calculation. 

One attractive feature of the cell model is that in the 
salt-free case the nonlinear PB equation can be solved 
exactly in this geometry |^. It also displays in a 
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clear fashion the effect of partial counterion condensa- 
tion 0, 111 • While a charged sphere loses all its counte- 
rions upon dilution, a charged plane keeps all of them. 
For a charged cylinder the fraction of ions, which upon 
dilution remain in the vicinity of the macroion, can be 
anywhere between and 100%. Addition of salt increases 
the screening of the charged rod. If the salt content is 
large enough, the electric field will have decayed to zero 
before the cell radius is reached. It is then permissible to 
extend the latter to infinity. This is the approach that we 
use for our integral equation calculations. 

The PB theory shows its limitations for systems where 
ion correlations become important 0. It cannot pre- 
dict the attractive forces which are seen experimentally 
in DNA solutions ^3 and which have been also found 
in various simulation 0, 0, 0| and integral equation 
[l^ nil nif studies in the presence of multivalent coun- 
terions. 

In this paper we present studies for a rod immersed in 
a salt solution which PB theory fails to describe correctly, 
namely, the possibility of overcharging a single rod, re- 
sulting in an effective charge reversal, and the appearance 
of a non-monotonic zeta-potential |l7lll8lll9l |. This effect 
is most relevant in electrophoresis experiments The 
overcharging of a macroion modifies in a non-linear fash- 
ion the electrophoretic mobility as a function of the poly- 
electrolyte charge and salt concentration. In the classical 
electrophoresis theory of Wiersema, O'Brien and White 
|2ll |22|. which is a linear theory based on the Poisson- 
Boltzmann (PB) description of the EDL around the poly- 
electrolyte, this effect is not included. 

The EDL obtained through the well established inte- 
gral equation approach does predict overcharging [r^ll8| . 
The overcharging effect has been observed in Monte Carlo 
(MC) simulations for planar [i^, cylindrical and 
spherical |2rI . geometries. Moreover, several years 
ago, Gonzales- Tovar et a/.^3| predicted that for some 
macroion charge and electrolyte concentration the over- 
charging produces an electrophoresis mobility reversal. 
This mobility reversal has been experimentally observed 
[iri li^ IsOl ■ This effect has been taken into account 
in the electrophoresis theory by Lozada-Cassou et al. . 
An excellent qualitative agreement with experimental re- 
sults ^0^28^,1^ is found. This overcharging effect and the 
predicted mobility reversal have been recently addressed 
by Deserno et a/.||j| and Shklovskii |33.l3^. 

In Ref. ion distributions are studied via the inte- 
grated charge distribution function, P(r), and overcharg- 
ing is addressed based on MD simulations. In the present 
paper we extend this research by approaching the over- 
charging effect through a well established integral equa- 
tion theory [itI as well as more detailed molecular 
dynamics (MD) simulations. We compute P{r), the con- 
centration profiles gi{r), the mean electrostatic potential 
profile ipir), and the ^-potential. A comparison of the in- 
tegral equation results with those from MD is presented. 




FIG. 1: Realization of the cell model. A rod of length Lb 
placed parallel to an edge of a cube of side length Lb yields 
an infinite square array of infinitely long rods upon periodic 
replication of the cubic box. The dark particles are the coun- 
terions or positive salt ions and the bright particles are the 
negative salt ions. This particular snapshot contains 36 diva- 
lent counterions, 220 positive and 220 negative divalent salt 
ions. 



SIMULATIONAL DETAILS OF THE MODEL 
SYSTEM 

Generating a cell-geometry 

Compared to the spherical cell model, the cylindrical 
one presents one additional but crucial complication: the 
charged rod is infinitely long. Several methods have been 
proposed in the literature to handle this problem. They 
essentially all use as a unit cell a hexagonal prism with 
a certain height. This approximates the cylindrical cell 
by the "most round" space-filling object. In this publi- 
cation we take a computationally even simpler approach 
by using a cubic unit box of side-length Lb and placing 
the DNA parallel to one of the edges (see Fig. ^ . The 
justification for this approach is as follows: In a high- 
salt environment the charged rod will be screened within 
a distance much smaller than the size of the box. The 
potential will therefore become constant long before the 
deviations from a cylindrical cell will be felt. It is thus 
largely irrelevant, whether these deviations are of hexag- 
onal or of cubic type. When comparing with the cylindri- 
cal cell model, a radius R — Li^/y/n is appropriate, since 
it yields the same area per rod. 

The main advantage of this approach is that such a 
system can be treated with the plain cubic Ewald sum 
or one of its mesh-upgrades 

m El 

This per- 
mits a very efficient way for computing the long-range 
electrostatic interactions. 
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Interaction potentials and DNA mapping 

A specification of two interaction potentials is neces- 
sary to describe the model system: (i) an excluded vol- 
ume interaction preventing two particles from occupy- 
ing the same position in space and (ii) the long-range 
Coulomb potential. 

For the excluded volume interaction between the ions 
we use the following potential: 
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Without the cutoff this would be the common Lennard- 
Jones potential describing particles with a diameter a and 
an attractive potential minimum of e. In order to achieve 
pure repulsion, we cut the potential at the minimum and 
shift it up such that it smoothly goes to zero there. As a 
consequence, e becomes largely irrelevant. Since two ions 
a distance a apart have the (repulsive) interaction energy 
e, we make the natural choice e = k^T. 

In order to give the rod a larger radius than the ions, we 
employ an ion-rod potential similar to Eq. in which r 
is replaced by r — rg. This shifts the hard core a distance 
Tg towards larger radii and gives a distance of closest 
approach of rg + a. Since the ions have diameter a this 
corresponds to a rod radius tq of tq = rs + a/2. Of course, 
in this case r is a cylindrical and not a spherical radial 
coordinate. 

The electrostatic interaction energy between two ions 
with charge z^eo and zjCo, is given by 



ZiZj£B 



with 



(2) 

where Zi is the valence of the i species and bq the ele- 
mental charge. By definition the Bjerrum length £b is 
the distance at which two unit charges have interac- 
tion energy /cbT'- With the vacuum dielectric constant 
Eq = 8.85 X lO-i^cV-^N-i and the dielectric coeffi- 
cient e-c = 78.5 applying to water at room temperature, 
one gets a Bjerrum length of approximately 7.14 A. The 
line charge density A of the rod is modeled by placing 
unit charges along the rod axis at the distance b = eo/X. 
The number of charges along the rod per Bjerrum length 
is an important dimensionless measure of the line charge 
density and is often referred to as the Manning parame- 
ter: 
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Its relevance lies in the fact that for ^ > 1 the phe- 
nomenon of counterion "condensation" is observed 

Within the periodic boundary conditions employed 
during the simulations, the presence of such long-range 



< r < rcut 
''cut < r. 



interactions poses both mathematical and technical diffi- 
culties. We use an efficient FFT accelerated Ewald sum, 
the P^M algorithm, which scales almost linearly with the 
number of charges [3^ |3a, |33 ■ 

The final step consists of explicitly mapping the pa- 
rameters to a DNA system in aqueous solution. This af- 
fects ion diameter a, rod radius tq, Bjerrum length £b and 
line charge density A. Our choice is presented in Tabic 

We performed Molecular-Dynamics simulations of the 
rod-systems using a Langevin thermostat in order to 
_d^\\f%the system into the canonical state [s^. The real- 
"space and Fourier-space part of the electrostatic energy 
was used to check for equilibration. The presented ob- 
servables originate from averaging over roughly 1500 in- 
dependent configurations. 



HNC/MSA THEORY 

The integral equation formalism is a well established 
statistical mechanical approach that has been shown 
to be reliable when it is applied to simple microscopic 
models of inhomogeneous fluids j33|, l40|] . A particularly 
successful integral equation theory for inhomogeneous, 
charged fluids is the so-called hypernetted-chain/mean 
spherical approximation (HNC/MSA) [H |41|, ,4^ 
The HNC/MSA theory for a model charged rod, im- 
mersed in an electrolyte, has been derived in the past 
[H nil. In this paper, we solve the HNC/MSA equation 
for an infinitely long, hard, charged cylinder of radius tq, 
with uniform line charge density A. The rod is immersed 
in a two-component restricted primitive model electrolyte 
(RPM), i.e., a fluid of charged hard spheres of diameter 
a with a centered point charge ZiEQ. 

The fluid electroneutrality condition is 



E 



iPm = 0, 



(4) 



where pm is the bulk concentration of species m. In order 
to satisfy the necessary condition of zero electrical field 
at infinity, the rod charge is compensated by the charge 
induced in the fluid. A': 
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-A, with pei(r) = eo ^ z„p„(r). 

m— 1 

(5) 

The local concentration profile of the species m is denoted 
as Pm{f)- The solvent is taken as a dielectric continuum 
of dielectric constant Er. For simplicity, rod and solvent 
are assumed to have the same dielectric constant, in order 
to avoid complications with dielectric boundaries. 

It has long been recognized in physics that particles 
and fields are equivalent in the sense that both are defined 
through their interaction potentials. This simple fact has 
been applied in the past to derive, in a straightforward 
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manner, inhomogeneous integral equation theories from 
the Ornstein-Zernike (OZ) equation for homogeneous flu- 
ids 0)01 • The homogeneous three-component OZ equa- 
tion is 



hij {r2i) 



,(1-21) 
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(r23)cmi(ri3) dv-i 



(6) 

where hij{v2i) and Cij(r2i) are the total and direct corre- 
lation functions, between particle 2 located at Y2 and par- 
ticle 1 at ri of species i and j, respectively, T21 = r2 — ri. 

In consequence, one can consider a particle of a fluid as 
a source of an external field or an external field as a par- 
ticle of the fluid. Applying this simple idea to Eqn. 
one can think of one of the species, say a, as made of 
infinite cylindrical rods and the remaining two species as 
ions. In a's infinite dilution limit the cylinders are un- 
correlated. Hence, letting pa — » 0, the Ornstein-Zernike 
equation for an ionic solution next to a charged cylinder 
reads 



2 . 

(r23)Cmj(ri3) 

m=l 

(7) 

In the past, several approximations (or "closures") for 
the direct correlation function have been suggested. Two 
of them are: 

ln(5ij(r2i)) = - PVij {r2i) + hij {r2i) - Cij {r2i) (8) 
Cy(r2i) = -/3Vy(r2i). (9) 

Eqn. (jHJ and Eqn. © arc known as the hypernetted chain 
(HNC) equation and the mean spherical approximation 
(MSA), respectively; Vij{r2i) is the direct interaction po- 
tential between species i and j, and /3 = l/fceT. 

If the HNC closure is used for the direct correlation 
function between the rod particle and the j species, equa- 
tion 121) becomes 



ffQj(i"2i) = exp{-/3V„j(r2i) 



(10) 
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In this scheme, the two-particles correlation functions 
haj{Y2i) and Caj{Y2i) correspond to the one-particle to- 
tal and direct inhomogeneous correlation functions hj (r) 
and Cj (r) for species j of a fluid under the influence of an 
external field produced by a rod particle. The local con- 
centration for the j species is given by Pj{v) = pjgj(r), 
where gj{r) = hj{r) + 1 is called the reduced concentra- 
tion profile. Therefore, the charge concentration profile 
in Eqn. © is given by 
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For this model the direct interaction potential between 
the rod and the j species of the fluid, Vj{r), can be sep- 
arated into two parts: the hard sphere-hard rod term 
Vj*{r) and the electrostatic interaction potential Vf^{r). 
The first takes into account the fact that ions cannot 
penetrate or deform the cylinder 

Too : < r <ro+ a/2 

\ : r>ro+a/2 
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The second can be found from Gauss' law and is given 

by 



-/3y/(r) =2z,ein(r) (r > ro). 



(13) 



In the MSA closure the homogeneous direct correlation 
function for a RPM electrolyte has an analytical expres- 
sion. This hmction can be written as 
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In the first term appears the direct electric interaction 
potential between the species of the fluid, which is given 
by Eqn. |(5J, the second term is an electrical short range 
function and the third is the direct correlation function 
for a hard sphere fluid. When the MSA closure is used 
in the integral of Eqn. one obtains the HNC/MSA 

equation for an ionic fluid next to a charged rod. Taking 
advantage of the cylindrical geometry and the fact that 
the direct correlation function between the ions depends 
only on their relative distance s = |ri — r.^l, Eqn. I|ll|l 
can after some algebra be written as [13, llal 



9]{r) = exp {-p{zjeoip{r) + Jj{r)) } , 
where ^{r) is the mean electrostatic potential 
- (3eoi]{r) = 2^1n(r) 



(15) 
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The Jj{r) terms are integrals of the short range terms of 
the direct correlation function 
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and 



L{r,y) = 
K{r,y) = 




l'ro+a/2 

A{r) ^ - K{r,y)dy 
Jq 

where = z'^ + r"^ + y'^ — 2ry cos 0. In the Uniit a 
of point ions, c^''(s) and c^''(s) ^ 0. Thus Jj (r) ^ [13 
and Eqn. p5|) becomes 

9m{r) = exp {-/Jz^eoV'lf') } , (19) 

where ip{r) is given by Eqn. p7|) . which is the solution 
of the weh known PB differential equation for point ions 
around a charged cylindrical electrode, i.e., 

l d^p{r) ^ d ip{r) ^ — eo_ ZmPmejip {~(3z„ieoil^{r) }. 
r ar dr^ ene,. -"^ — ' 

m— 1 

(20) 

Therefore, the integral equation version of the Poisson- 
Boltzmann differential equation is [T^ 

.gj(r) =exp{2z,ein(rX21) 
+ 2z,7r^B / Pcd(y) ln((r' +y^ + - y^\)/2)ydy 1(^2) 

Jro+a/2 J 

which is also the HNC/MSA equation in the point ions 
limit. The differences between these two theories are due 
to the ionic size correlations which are partially taken 
into account in HNC /MSA theory but which are ignored 
in the PB equation. In the past, it has been shown that 
the size effects become important for strong field inter- 
action, poly-valent ions and high salt concentration |l8l| . 

The HNC /MSA and PB integral equations are numeri- 
cally solved with efficient finite element methods IQn^^- 
The solution of Eqn. HI5|I takes one minute in a R12000 
processor of an SGI machine. In our theoretical calcu- 
lation we use the same input parameters as in our MD 
calculations. The salt concentration is obtained through 

Ps = N,/Ll (23) 

where Ng is the number of salt molecules used in the sim- 
ulation. Nevertheless, we have to point out the differences 
in the short range interaction used in both models. While 
in HNC/MSA we use hard particles, the MD for technical 
reasons uses the potential from eq.^ which has some sur- 
face softness. This results in an effective excluded volume, 
whose actual value depends on the interaction strength. 
At an interaction energy of /cbT the particle diameter is 
a. However, if the particles attract each other strongly, 
they can come closer to each other. 



parameter 


symbol 


value 


value in LJ units 


ion diameter 


(X 


4.25 A 


d 


ion valence 


z 


2 


2 




To 


7.86 A 


1.85 a 


line charge density (DNA) 


A 


-eo/l.7 A 


— 2.5 Co / 0, 


Bjerrum length, (water) 


£b 


7.14 A 


1.68 a 


Manning parameter 


C 


4.2 


4.2 


box size 


ib 


122.4 A 


28.8 a 


corresponding cell radius 


R 


69.1 A 


16.2 a 


temperature 


T 


298 K 





TABLE I: System parameters for the DNA simulations. 



RESULTS 

We will discuss the results in terms of the mean elec- 
trostatic potential, ipir), the local concentration profiles, 
gmiT)^ and the integrated charge distribution function 

= m / ^r2nfp,i{f). (24) 

I'^l Jra 

Figure H shows the HNC/MSA and MD results for P{r) 
and tl^{r), for five rod-systems, for which the line charge 
density, and thereby the Manning parameter, has been 
successively increased. One system has a smaller Man- 
ning parameter than DNA and three a larger. All sys- 
tems have a Debye screening length which is much smaller 
than the Manning radius of the corresponding salt free 
system. According to the observations of Ref. this 
means that the concept of Manning "condensation" is no 
longer meaningful. The rod charge is compensated by salt 
screening. According to simple PB theory, no remarkable 
features are to be expected in the bulk solution. Con- 
trary to that, however, in Fig.|3^ P{f) shows overcharg- 
ing, i.e., the rod charge is over-compensated at a certain 
distance from the rod. For a given salt concentration the 
degree of overcharging is increased as the charge density 
increases. No overcharging saturation point is observed, 
at least for physical charge densities. While in Fig. |2K 
overcharging is present even for a very low charge den- 
sity, such as a = COOSC/m^ our HNC/MSA calculations 
show overcharging for even lower rod charge densities. 
We will come back to this point later. In Fig. Oi, the 
maximum of overcharging is closer to the rod's suface as 
a increases. This effect has implication for the location 
of the so-called C— potential, in electrophoresis experi- 
ments [13 ■ As pointed out before 0,1131, overcharging 
implies a change of direction in the local electrical field. 
Our MD and HNC/MSA results are consistent with this 
fact, as shown in Fig. ^3p, where a minimum in '4>(r) is 
observed. The value of r at the minimum of V'('')j i-S-, 
where the electrical field is zero, concides with that at 
which P{r) = 1, as it should. 

Integral equation theories give much more reliable re- 
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C Clpb [mV] Cpb [mV] (hnc [mV] Cmd [mV] 


















1.05 


10.52 


10.28 


6.262 


9.116 


2.1 


21.04 


19.48 


10.99 


14.62 


4.2 


42.09 


33.32 


14.27 


16.25 


6.3 


63.13 


42.96 


11.08 


13.80 


8.4 


84.17 


50.16 


3.771 


9.263 


10.5 


105.2 


56.03 


-6.219 


3.675 



TABLE 11: C-potential of a DNA-sized rod (see Tab. |TJ im- 
mersed into a 0.49 M electrolyte of 2:2 salt as a function of 
its Manning parameter ^. The four predictions are from lin- 
ear PB theory, PB theory, hypernetted-chain theory and MD 
simulation. The error in the latter is estimated to be of the 
order of 2%. Figure [^visualizes the data. 



suits than Poisson-Boltzmann at high salt concentra- 
tions, poly-valent ions and/or for fluids under the in- 
fluence of strong external fields. This is because they 
partially take into account correlations due to the ionic 
size, which PB theory ignores. It can be seen that 
the HNC/MSA theory indeed correctly predicts the 
occurrence of overcharging, although it overestimates 
its amount. The qualitative agreement between the 
HNC/MSA and MD results is excellent. The quantita- 
tive agreement is fair and particularly good for the DNA 
calculation. In the next section we will discuss the dis- 
agreement seen between the HNC /MSA and MD results, 
for high and low cylinder charge densities. 

For a given constant value of the Manning parameter, 
cylindrical macroions with different diameters will pro- 
duce different EDL structures. In Figure El HNC/MSA 
results for rods with the same value of ^ but different 
diameters show different charge distributions P{r). This 
illustrates the simple fact that the Manning parameter f 
alone is not sufficient to specify the ionic structure around 
cylindrical macroions. Since the ionic structure depends 
directly on the field strength, it is often better to use 
the cylinder's surface charge density a, which is a direct 
measure of the electrical field that penetrates into in the 
solution. Note, however, that ^ and a are simply related 
by ^ = 27rro(T€B/eo. 

In Figs. Hi andEb we show HNC/MSA and MD re- 
duced concentration profiles for a — 0.190 C/m^ = 
4.2) and a = 0.49 C/m^ ( ^ = 10.5), respectively, whereas 
all the other parameters are the same as in Fig. |21 The 
value of r at the maximum of P(r) matches with that 
where g+{r) = 9-{r), at least for symetric electrolytes. 
Therefore the later oscillation in the local concentration 
profiles, show the attraction (repulsion) to colons (coun- 
terions), due to a change in direction of the effective 
electrical field. On the other hand, concentration profile 
oscillations are a consequence of the ionic size correla- 
tions. Hence, overcharging results of electrostatic attrac- 
tion and size correlations. 
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FIG. 2: (a)Charge distribution function P(r) and (b) mean 
electrostatic potential i/'(r) for the parameters from Tab. 
The electrolyte is a 2:2 salt of concentration 0.49 M. The dis- 
tance r is measured with respect to the cylinder surface in 
units of the ionic radius. The solid and dashed lines corres- 
pond to the HNC/MSA and MD calculations, respectively. 
The different curves are for different values of the Manning 
parameter, ^ = 2.1, 4.2, 6.3, 8.4, 10.5, which correspond to sur- 
face charge densities of cr = 0.095,0.190,0.286,0.381,0.476 
C/m^, respectively. In P(r), the surface charge density in- 
creases from bottom to top, whereas in the mean electrostatic 
potential it increases from top to bottom. The value ^ = 4.2 
(cr = 0.190 C/m^) corresponds to DNA. 



The first maximum seen in the colon reduced concen- 
tration profile implies a colon concentration above its 
bulk value and indicates that the local electric field is 
attractive to colons. It is also observed that HNC/MSA 
overestimates the contact value of the distribution func- 
tion with respect to molecular dynamics predictions. The 
overestimation of HNC /MSA can be associated with two 
facts: (i) the difference in excluded volume used in both 
models; and (ii) HNC/MSA theory does not take into 
account all the size and charge correlations. 

Another important phenomenon is indicated by the be- 
havior of the mean electrostatic potential in Fig. It 
concerns the value of this potential at the distance of 



7 



1.2 




r[a/2] 



FIG. 3: Charge distribution function from HNC/MSA, for a 
RPM electrolyte next to a charged cylinder. The distance r 
is measured with respect to the cylinder surface in units of 
the ionic radius. The electrolyte is a 2:2, 0.49 M salt. The re- 
sults correspond to different cylinder radii, whereas the Man- 
ning parameter has been held constant at ^ = 4.2. From top 
to bottom the cylinder radii are ro = 3.86, 7.86, 110, 1200 
A, which correspond to the surface charge densities a = 
0.388,0.190,10-^,10"^ C/m^, respectively. The curves for 
ro = 110 A and ro = 1200A are indistinguishable on this scale. 

closest approach between ions and the rod, 

C^^{ro + a/2), (25) 

which can be identified with the zeta-potential of elec- 
trophoresis and which is of importance in the compu- 
tation of electrophoretic mobihties [Ulil. In the PB 
theory for a charged rod immersed into an ionic solu- 
tion, and its linearized version, ^ depends monotonically 
on ^ [13] ■ However, Fig. |21 shows that, in the presence of 
overcharging, ■!/)(r) can become negative. This calls the 
monotonic relation between ^ and ^ into question. Ta- 
ble ^1 summarizes various predictions for the C-potential 
as a function of Manning parameter together with results 
from the molecular dynamics simulations. A graphical il- 
lustration is given in Fig. (S) 

Indeed, C, is found to first increase with ^, but from 
a certain value on it decreases and would finally even 
become negative. Note that this would reverse the drift 
direction in electrophoresis measurements, as first pre- 
dicted by Gonzales- Tovar et aL0| and recently demon- 
strated, for spherical macroions, by Lozada-Cassou et 
al.^^. While nonlinear and linearized PB theory coin- 
cide with the data and with each other for small Man- 
ning parameter, they completely fail to predict the back- 
bending, which already sets in at comparatively small 
values of The HNC/MSA theory captures this effect, 
but it underestimates the value of the potential. How- 
ever, comparisons for the C-potential, as a function of cr, 
between HNC/MSA and MC simulations for planar and 
spherical hard, charged macroions, immersed into a RPM 




r[a/2] 

FIG. 4: Reduced concentration profiles gi{r) for an electrolyte 
next to a charged cylinder. The electrolyte is a 2:2, 0.49 M 
salt. The distance r is measured with respect to the cylinder 
surface in units of the ionic radius. The system parameters 
are given in Tab. |I] Fig. 2t shows the results corresponding 
to ^ = 4.2 (cr = 0.190 C/m^) and Fig. gja the resuhs for 
^ = 10.5 (cr = 0.49 C/m^). The dots and crosses are the MD 
distribution functions for counterions and colons, respectively. 
The sohd lines are the resuhs from HNC/MSA. 

electrolyte, have been made in the past URESES- For 
these two geometries, for 1:1 electrolytes and for 2:2, at 
around 0.5 M, as in our Fig.El the HNC/MSA-MC agree- 
ment is excellent. The disagreement seen in Fig.Elcan be 
related to the different short-range potentials used in our 
MD and HNC/MSA calculations, which gives for the MD 
not a sharply defined ion-cylinder contact. This is impor- 
tant, since the ^-potential is evaluated exactly there. For 
planar and spherical macroions, PB C vs. charge curves 
are also above the HNC/MSA and MC curves, which is 
consistent with our present observations. 

For a given rod charge and diameter, for any electrolyte 
solution parameters, a higher C-POtential always implies a 
lower counterions adsorption A soft-repulsive short- 
range potential, as in our MD calculations, allows the 
counterions to get closer to the rod. This produces a 
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FIG. 5: Zeta potential C = i/'(ro + a/2) as a function of surface 
charge density o", for a cylinder immersed in a 0.49 M, 2:2 
electrolyte. The model parameters are given in Tab. ^ The 
dotted line is the prediction of PB theory, the dots are the 
MD results and the solid line is the results from HNC/MSA 
calculations. 



higher ^-potential in the MD simulation when compared 
to the HNC/MSA results. In the case of the PB theory, 
the neglect of the ionic-size correlations allows a higher 
density of counterions next to the rod, which in turn can 
never lead to overcharging. 

In all cases, the HNC/MSA qualitative agreement with 
MD is very good. The success of HNC/MSA theory in- 
dicates that local ion-size correlations are responsible for 
both phenomena. As it is concluded by Lozada-Cassou 
and Jimenez- Angeles in ref. |4^ . overcharging increases 
with an increase of the system excluded volume: i.e., 
higher ionic size and/or concentration. If size correlations 
are not considered (as in PB theory) or are negligible, 
overcharging does not occur, even for divalent ions. On 
the other hand, overcharging can occur in a high excluded 
volume system even for monovalent ions. To illustrate 
this fact we have simulated two systems, for which the 
distance of closest approach to the charged rod is in both 
cases 9.98 A, hence DNA-like. The first rod system has a 
charge parameter of ^ = 4.2 (a = 0.240 C/m^) and is im- 
mersed in a 1.0 Molar 1:1 electrolyte solution with ions of 
diameter a = 7.4 A. The second rod system has an even 
higher charge parameter of ^ = 10.5 ( cr = 0.40 C/m^) 
and is immersed in a 0.5 Molar 2:2 electrolyte, but, this 
time, the diameter is only a = 1.0 A. The results are 
shown in Fig. |S1 where the charge distribution function 
is plotted. Contrary to the general believe, the system 
with the large monovalent ions shows overcharging, both 
in the HNC/MSA and the MD results, and the location of 
the peaks coincide, whereas we see small differences in the 
height. Further, strikingly, the system with small divalent 
ions shows no sign of overcharging, although the charge 



FIG. 6: Charge distribution function from HNC/MSA (solid 
lines) and MD (dotted lines), for a RPM electrolyte next to a 
charged cylinder to which the ions have a distance of closest 
approach of 9.98A. The distance r is measured with respect 
to the cylinder surface in units of the ionic radius. In this 
plot a comparison between a low and a high excluded volume 
systems is presented. The system where overcharging occurs 
corresponds to a monovalent electrolyte with p = 1.0 M, ^ = 
4.2 (a = 0.240 C/m^), and large ions (a = 7.4A). The system 
where overcharging is not observed corresponds to a divalent 
electrolyte with p = 0.5M, ^ = 10.5 {a = 0.40 C/m^), and 
small ions (a = I.OA). 

parameter for this system is even higher. The shoulder 
in the distribution of the MD simulation indicates that 
the second ionic layer is almost neutral. A possible ex- 
planation is the following: The small ions have an energy 
of attraction of almost SO/cb^, hence the salt is actually 
present in the form of little overall neutral salt clusters. 
This is supported by visual inspection of the MD con- 
figurations. These mostly neutral clusters are polarizable 
and hence tend to accumulate in the large field gradient 
surrounding the rod. The free counterions left can come 
only from the rods charge and can in no way lead to an 
overcharging. A key prerequisite for this scenario is that 
the interaction energy between ions exceeds the interac- 
tion between ions and the rod. Also consistent is, that 
the HNC/MSA calculation as well shows the vanishing 
of the overcharging. 

Finally, one might ask the question of how the ef- 
fect of overcharging depends on salt concentration. Ob- 
viously there can not be any overcharging in the ab- 
sence of salt, since then even a complete "condensation" 
of all ions would merely neutralize the rod. However, 
it is not clear from the beginning whether an addition 
of just some salt will immediately lead to overcharging. 
Figure □ shows HNC/MSA and MD distribution func- 
tions P{r) and mean electrostatic potentials tp{r) for a 
DNA-sized and charged rod, immersed into a 2:2 elec- 
trolyte, for varying salt concentrations. The phenomenon 
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FIG. 7: Charge distribution function P(r) (left) and mean 
electrostatic potential i/'(r) (right) results, for the system from 
Tab. ^ with Manning parameter ^ = 4.2. The distance r is 
measured with respect to the cylinder surface in units of the 
ionic radius. The line styles have the same meaning as in 
Fig. 121 The curves correspond to different salt concentration, 
Ps = 0.12,0.24,0.34,0.49,0.68 M. In the charge distribution 
function the salt content increases from bottom to top, in the 
mean electrostatic potential it increases from top to bottom. 

of overcharging is indeed observed, but only at sufRciently 
high salt concentration. These results are consistent with 
the fact that higher electrolyte concentration increases 
the excluded volume, thus increases overcharging |48| . 
Therefore, as proposed in Ref. 48], higher excluded vol- 
ume implies higher overcharging. This may be also illus- 
trated by defining ri to be the radius at which overcharg- 
ing sets in, i.e., 

ri = min {r : P{r) — l}. (26) 

Fig. IHl illustrates the measured ri together with a 
hypernetted-chain prediction as a function of salt concen- 
tration ps- Clearly, ri must increase with decreasing ps, 
since overcharging is reduced and consequently the size 
of the charge-compensating layer must increase. In fact, 
hypernetted-chain theory predicts that ri(ps) diverges at 
some finite density ~ 0.18mol/l corresponding to a 




0.1 0.3 0.5 0.7 

P/M7 

FIG. 8: Radius ri from Eqn. I|26^ at which overcharging sets 
in, for a DNA-sized and -charged cylinder (ro = 7.86A and 
a = 0.190 C/m^), as a function of salt concentration ps. The 
distance ri is measured with respect to the cylinder surface in 
units of the ionic radius. The dots are the results of molecular 
dynamics simulations while the solid line is the prediction of 
HNC/MSA. 

salt Debye length of 3.59 A, or roughly 200 salt molecules 
within the simulation box. For the lower concentrations, 
the simulated ri values lie below the hypernetted-chain 
prediction, but this is not a feature generally to be ex- 
pected. In Fig. [71 it can be seen that the stronger over- 
charging in hypernetted-chain theory should normally 
lead to a value of ri smaller than in the simulation. But 
with decreasing density the finite radius of the simulation 
cell becomes relevant. Indeed, zero-salt distribution func- 
tions reach the value 1 at the cell boundary, and not at 
infinity. A reduction of the amount of added salt at fixed 
cell radius must necessarily lead to values of ri smaller 
than the hypernetted chain prediction for the bulk, since 
within a finite cell ri cannot diverge. On the other hand, 
for technical reasons, the MD cell radius can not be in- 
creased indefinitely. In any case, our MD results have not 
been able to unambiguously detect overcharging at den- 
sities equal to or lower than 0.2 M, which is in excellent 
agreement with the HNC/MSA prediction for the mini- 
mum amount of salt needed to induce overcharging. From 
Figs. lO and (jEJ it is clear that excluded volume effects 
determine the occurrence of overcharging, as proposed by 
some of us ;48j . 



CONCLUSION 

Theoretical and numerical studies of stiff linear poly- 
electrolytes, immersed into a 2:2 RPM electrolyte have 
been presented. The qualitative agreement between the 
MD and HNC/MSA resuhs is excellent. For the par- 
ticular case of DNA parameters, there is also a very 
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good quantitative agreement. We argued that the bet- 
ter agreement in the planar and spherical case reported 
in Refs. [2^ 0, 0| is most probably due to an imple- 
mentation of the short range interactions that is identi- 
cal to the theoretical model and not merely very close. 
At low salt concentrations the difference between a fi- 
nite cell used in simulations and the R = 00 system 
used in the HNC/MSA calculations becomes important. 
Also, the macroion and ions electrical fields are much 
less screened and, hence, the charge correlations become 
more important [2| • Both these points are responsible for 
the disagreement between theory and simulations at low 
salt. While the former can in principle be resolved by in- 
creasing the size of the simulation cell, the latter can be 
addressed by the three point extension HNC/MSA theory 
^li^l, which gives a better account of these correlations 

In the past the HNC/MSA theory for an electrolyte 
next to a charged rod has been derived jl3| and applied 
(isf to simple DNA models. In this paper we have shown 
this theory to be qualitatively correct. In particular, for 
DNA parameters the quantitative agreement is excellent. 
To the best of our knowledge, in Gonzales- Tovar et oi.fl^ 
overcharging has first been predicted theoretically. It was 
point out that this should entail an electrophoresis mo- 
bility reversal. In the PB theory the ^-potential as a func- 
tion of the cylinder charge density is a monotonic func- 
tion, while in the HNC/MSA theory it is not. In this 
paper we have shown the HNC/MSA prediction to be 
in agreement with MD calculations. This non-monotonic 
behavior of the ^-potential has important implications in 
electrophoresis calculations. The standard electrophore- 
sis theory 0, is based on the PB prediction for ( 
and, hence, important differences should be found if the 
HNC /MSA theory is applied to the electrophoresis prob- 
lem. 

Recently, Lozada-Cassou et a/.jsj] extended the elec- 
trophoresis theory to include ionic size effects, through 
the HNC/MSA theory. It was applied to spherical 
macroions. A mobility reversal and a very non-linear be- 
havior of the mobility, as a function of the ^-potential, 
were predicted. These predictions were found to be in 
agreement with experimental results. Contrary to the 
standard electrophoresis theory, the mobility was shown 
to be non-universal. This is due, precisely, to the non- 
monotonic behavior of the ^-potential. The HNC/MSA 
and MC simulations for spherical and planar macroions 
do not show a maximum in C(cr) as pronounced as that 
reported here for the cylindrical geometry. Therefore, we 
believe that our HNC/MSA and MD results are particu- 
larly relevant. 

In view of the qualitative agreement between 
HNC/MSA and MD and the quantitative and qualitative 
disagreement of these with Poisson-Boltzmann theory, we 
conclude that for concentrated and divalent electrolytes, 
influenced by charged polyelectrolytes, the size correla- 



tions must be taken properly into account to describe 
the cylindrical double layer. Our MD results clearly in- 
dicate that overcharging strongly depends on the system 
excluded volume, as proposed by some of us While 
maximum overcharging does not seem to have a limit, 
the minimum conditions to have overcharging depend on 
many different ways in which ionic charge, size and con- 
centration, and surface charge density participate in a 
system, as we showed in Figs. ® and (0 
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